Imputation methods for missing failure times in recurrent-event survival analysis: Application to suicide attempts in the transgender population

Suicide risk among transgender populations is an important public health issue. In a project evaluating association between gender affirmation and suicide attempts in the US Transgender Survey, we evaluated the relationship between gender affirmation and risk for suicide attempts. One of the challenges is that the age at suicide attempts was only collected for the first and last attempt. The initial zero-inflated negative binomial model enabled us to evaluate the association between gender affirmation and number of suicide attempts per 5 years adjusting for other covariates. However, ignoring missing failure times of recurrent events may have caused bias and loss of efficiency. In this paper, we use a recurrent-event survival analysis incorporating time-varying covariates with three approaches to impute the age at suicide attempt, estimates from three imputation approaches are similar. We were able to confirm the findings from the initial model and identify additional associations that were not detected in the initial analysis. Findings suggest the need to consider additional analytical approaches in settings with high data missingness by design. Research to validate and compare measures that ask first and last attempt to those which enumerate all attempts in this population will be important for future surveys.


Introduction
Suicide risk among transgender youth and adult populations is an important public health issue [1][2][3]. Transgender individuals (those who have a gender identity that differs from their sex assigned at birth) have a two-to four-fold increased risk of suicide attempt compared to cisgender (non-transgender) individuals [4][5][6]. Identifying culturally-specific factors associated with suicide risk in vulnerable populations is important to inform national prevention strategies [7]. Gender affirmation is one such transgender-specific factor to consider in suicide risk. Gender affirmation refers to changing one's gender presentation to affirm or reflect one's internal gender identity [8]. It includes social affirmation (socially presenting in one's identified gender identity) such as using a selected name and pronoun, and medical affirmation (interventions or treatments that physically alter the body) such as hormones or surgery [9]. In 2015, the US Transgender Survey [10] was conducted to learn about transgender experiences on a variety of health topics, including suicide risk. A question of interest is whether there is a relationship between a gender affirmation and risk for suicide attempts. One part of the survey asked participants about suicide attempts. These questions provide information on the date of the first attempt and the date of the last attempt as well the number of attempts in between these two dates. The date of the attempts between the first and last date were not collected. If the times of each attempt were collected, a recurrent survival analysis could be implemented.
Due to the missing information on timing of each suicide attempt, the initial analysis [11,12] used a weighted time-varying survival analysis to evaluate the association of gender affirmation with time to first suicide and a zero-inflated negative binomial model with a 5-year offset and over-dispersion parameter to assess the association of suicide attempt rate over lifetime with gender affirmation. Both models were stratified by age at survey completion and were adjusted for other factors such as race, education, and age at transgender awareness.
While these approaches identified important associations, ignoring missing failure times of recurrent events may have caused bias and loss of efficiency. In the survival analysis of the first attempt, age at the first suicide and gender affirmation (time-varying covariate) were analyzed while information about the total number of attempts and age at the last attempt was discarded. Whereas, in the analysis of total number of attempts using a zero-inflated negative binomial model, number of attempts offset by total years of follow-up and gender affirmation as static variables were analyzed, while information on timing of attempts and gender affirmation was discarded, causing lack of temporality between exposures and outcome. To overcome the shortcomings in both analyses, it is important to consider alternative approaches that would incorporate all the information on exposure and outcome simultaneously.
The missing suicide attempt dates can be viewed as a missing data problem. In this case the missing data provide information on the timing of the outcome. Missing data is a common challenge in research studies and a large body of literature exists addressing the statistical issues that arise [13][14][15]. It is well recognized that missing data results in bias and efficiency loss and several methods have been proposed to address this issue [16][17][18][19][20]. Multiple imputation is one of the approaches that is often used to address this problem. Briefly, multiple imputation involves generating multiple datasets by replacing the missing values in the sample data with plausible draws generated from the predictive distribution of models fit on the observed data. One of the assumptions to implement multiple imputation is that the data are missing completely at random (MCAR) or at least missing at random (MAR) [13]. The data are said to be MCAR if the probability of being missing is the same for all cases. The data are said to be MAR if missingness is systematically related to the observed data but not the unobserved [21]. Given that the missingness in this study is by design, it does not depend on the values would have been observed. Thus, it is reasonable to assume that these data are at least missing at random and therefore, conventional imputation methods can be utilized. For the subjects with more than two suicide attempts we do not have the age of attempt, but we do know this age fell between the age at the first and the last attempt. If we had the exact age of every suicide attempt, we could utilize recurrent survival analysis methods. In this paper, we implement three imputation approaches utilizing the age distribution of the suicide attempt and compare the results of recurrent survival analysis for the three approaches as well as to the results from the initial zero-inflated negative binomial model that did not involve missing data imputation. If the results are different this could result in some recommendations for modifying the study design for future surveys.

Study population
The 2015 U.S. Transgender Survey (USTS) is the largest survey targeted at the transgender adults residing in the United States. Details on survey methodology are published elsewhere [11]. "transgender" was defined broadly for this study to include any individuals who identify as transgender, trans, genderqueer, non-binary, and other identities on the transgender identity spectrum. The web-based survey had 32 sections and 324 possible questions that covered various topics. The study sample included 27,715 respondents for whom data were collected over a 34-day period in the summer of 2015, among which 24,343 reported wanting or having had surgery or hormones. There were 4561 participants excluded because a) age of transgender awareness or age at first attempt was missing or b) age of transgender awareness was after date of suicide attempt or date of gender affirmation, or last attempt. Among the remaining 19,782, there are 14,231 participants who identify as binary. Analyses were a priori restricted to binary-identified participants (i.e., trans men, trans women) because gender affirmation needs and practices differ for this group relative to non-binary individuals (i.e., genderqueer, genderfluid) [22].
Written informed consent was obtained for all subjects included in this study. The survey used informed consent with a waiver of signature to protect anonymity of survey subjects. The participants were given a study info sheet at the beginning of the survey and were asked if they agreed and consented to take the survey before proceeding into the survey. The survey was reviewed and approved by the UCLA North General IRB.

Important covariates
Various covariates are used in the analysis including gender affirmation, age at survey completion (current age), age at transgender awareness, race, ethnicity, assigned sex at birth, and education [11].
Three types of gender affirmation processes are of interest: social affirmation, hormones, and surgery. While the survey captures age at the first uptake of each gender affirmation process, it does not provide information on dosage and types of hormones, or number and types of surgical procedures received over gender time. Thus, it is difficult to study dose-response precisely. Based on age at the first use of each type of gender affirmation, we derived three time-varying binary variables indicating whether the participant had certain gender affirmation at a specific time (denoted as t). Age at survey completion was grouped into four categories: 18-24, 25-29, 30-39, and 40 or above. Age at transgender awareness (first recognizing oneself as being transgender) was also grouped into four categories: 10 or below, 10-14 (including 14), 14-18 (including 18), and above 18. We categorized these two age variables differently because awareness usually occurs at an early age. Race questions were drawn from the American Community Survey (ACS) and collapsed into 6 categories: Alaska Native/American Indian alone, Asian/Native Hawaiians (NH)/Pacific Islanders (PI), Biracial/Multiracial/Not listed, Black/African American alone, Latino/a/Hispanic alone, and White/Middle Eastern (ME)/ North African (NA) alone. Sex assigned at birth (male or female) was included as a binary variable. The highest level of education was recorded into 6 categories based on ACS: less than high school, high school graduate (including GED), some college (no degree), associate's degree, bachelor's degree, and graduate or professional degree.

Assessment of suicide attempts
The outcome of interest for our analysis is lifetime suicide attempts. Several variables were derived from the survey: 1) a binary indicator whether the participant had attempted suicide, 2) the number of attempts, 3) age at first attempt, and 4) age at last attempt. The total number of attempts was truncated at 26. If there were additional attempts between the first and last attempt, we do not know at what age they occurred. Consequently, recurrent-event survival analysis, which would otherwise be a natural choice for this type of data, becomes impossible to conduct directly. To illustrate the challenge in analyzing these data when age at suicide attempt is missing for some participants at some events, Fig 1 compares two participants having different numbers of attempts and ages. As shown, participant (a) had two attempts and no missing data because there are no attempts before the first and last attempts. Therefore, temporal relationship between attempts and gender affirmation is clear, making the assessment of the association between gender affirmation and suicide attempts relatively straightforward. By contrast, participant (b) had five suicide attempts and the age at the 2 nd , 3 rd , and 4 th attempt are unknown. The arrows show all possible timings of these events relative to exposures, which suggest different associations between exposures and outcomes.

Statistical analysis
In order to analyze all the attempts and evaluate their temporality relative to exposures, we impute missing age at attempt using three approaches. For each imputation approach the analysis consists of three steps: 1. Generating imputed data sets 2. Analyzing each of the imputed data sets using recurrent event survival analysis 3. Combining the estimates of interest using Rubin's rules [13].
The following sections provide details on each step.

Imputation methods
The three approaches differ on how much the imputation depends on the observed data. By comparing results from these approaches, we can then evaluate the impact of the assumptions. First, we assume age at attempt follows a uniform distribution irrelevant to the observed data. Next, we allow a slightly stronger dependence on the observed data by assuming the distribution of attempts over age is the same in the observed and unobserved data. Last, we assume covariates also affect the distribution of age at attempt, in the same manner in the observed and unobserved data. Details on assumptions and implementations are discussed below.
First, a simple random imputation (denoted as SRI) approach assumes attempts randomly occurred between age at first attempt and age at last attempt following a uniform distribution (approach 1). An integer was draw with replacement from a uniform distribution ranging from age at first attempt to age at last attempt for each attempt with missing age.
Secondly, a modified SRI approach assumes attempts randomly occurred between age at first attempt and age at last attempt with probabilities of the observed values of age at attempt used to generate the missing data (approach 2). There are three steps to this approach. First, the probabilities of making suicide attempts at all ages were obtained from the observed dataset with each participant contributing up to two ages. Second, for each unique pair of age at first attempt and age at last attempt, probabilities within the range were standardized to add up to 1. Finally, an integer was drawn with replacement from a distribution using the probabilities from step 2 for each attempt with missing age. For each approach, 10 imputation data sets were generated.
The two approaches above assumed missing completely at random (MCAR) within age stratum while the data may actually be missing at random (MAR). To accommodate this possibility, a multiple imputation (MI) approach (approach 3) using linear regression method for continuous variable with univariate missing pattern was also explored. Data were transposed so that each row was an attempt and year of attempt is the only variable that could be missing. Age at attempt was assumed to be normally distributed and the imputed values were rounded to the nearest integers. Variables used to predict missing ages included age at first attempt, age at last attempt, age at survey, integer representing the order of attempt (2 nd , 3 rd , etc.), age at transgender awareness, race/ethnicity, sex, and maximum level of education. Note that age at survey can be treated as a baseline characteristic and included in the regression model because all the participants were surveyed at the same year and age at survey is equivalent to year of birth, which was prior to age at transgender awareness. To impose participant-specific upper bound (age at last attempt) and lower bound (age at first attempt), after each imputation, values out of range were re-set to missing and the imputation process was repeated until all values were within range (this required <100 iterations). Attempts of each participant were sorted by age and a new order of attempt was generated to replace the old. We also imputed on the logtransformed scale due to observed right-skewness and compared it to the original to assess the sensitivity of results to normality assumption. Here we rounded age to be consistent with the other two approaches; however, results are similar if we do not round age. We used the monotone regression method implemented by the standard SAS procedure PROC MI with MONO-TONE REG. Details of this method was documented elsewhere [23,24].
We implemented all three approaches and examined how the modelling results differed.

Recurrent event survival analysis
After all ages of suicide attempts were imputed as above, the next step is to evaluate the association of time to suicide attempt and time-varying gender affirmation using recurrent event survival analysis. Modeling was carried out using a Cox proportional-hazards (PH) model and it is assumed that a) all the recurrent events are identical, i.e., the order of event is negligible and b) the PH assumption is satisfied for all variables. Unlike non-recurrent event models, this recurrent event model kept the participants in the risk set until the last attempt or censoring. The data were constructed so that each observation, or interval, corresponds to each recurrent event while gender affirmation remained unchanged. For participants with more than one interval, the different observations contributed by the same participants were treated as if they were independent. To adjust for the likely correlation among recurrent events on the same participant, robust variance estimates were used. All the analyses were conducted in SAS 9.4. The model can be written as: Gender affirmation, represented by three variables including social affirmation, hormones, and surgery, was considered time-dependent, as denoted by the X j (t) variables. All the other variables, including gender, age at transgender awareness, race/ethnicity, and maximum level of education, were considered time-independent, as denoted by the X i variables. Maximum education acted as a proxy for social economic status of the participant. Thus, it was preferred over timevarying level of education, which was confounded by lapse of time since transgender awareness.

Combining the estimates
The recurrent event survival analysis was conducted on each imputed dataset using the PROC PHREG procedure. Then parameter estimates and associated variance estimates were combined using the PROC MIANALYZE procedure to derive valid inferences for these parameters.
With m imputations, the combined point estimate of a parameter is the mean of the m individual estimates. By Rubin's rule, the combined variance is: Where � W is the within-imputation variance of the variance estimates, and B is the between-imputation variance of the variance estimates.
Simulation studies show that this approach has reasonable coverage, bias and mean squared error (S4 Table).

Results
The demographics are summarized in the initial report [11]. From this report, the study population is 84% white, 57% male assigned at birth and 31%, 18%, 20%, and 31% were 18-24, 25-29, 30-39, and > = 40 years of age, respectively, at survey completion. In addition, 35%, 25%, 19%, and 21% were < = 10, 11-14, 15-18 and >18 years of age, respectively, at awareness and 2%, 11%, 36%, 10%, 26%, and 15% of participants were in the following 6 education level categories: less than high school, high graduate, some college, associate's degree, bachelor's degree, and graduate/professional degree. Histograms (Fig 2) show the distributions of age at attempts in the 10 imputed datasets combined compared to the original data. Imputed datasets all have the distribution of age at suicide attempt similar to that in the originally observed data.
In the binary gender identity group, we found moderate associations of surgery and hormones on estimated hazard ratio (HR) as well as interactions with social affirmation that are differential by age. Table 1 shows the estimated unadjusted hazard ratio and the 95% confidence interval for each imputation approach. Table 2 and Fig 3 show the estimated adjusted HR and 95% confidence interval for each imputation approach. Surgery with social affirmation is consistently associated with decreased risk for suicide attempt across age groups, with small variation among imputation approaches ( Table 2,

PLOS ONE
(IRRs) from the initial analysis using zero-inflated negative binomial model are also included for comparison.
Simple random imputation assuming uniform distribution or using probabilities lead to similar point estimates and conclusions. The estimated HRs are consistent across imputation methods for most of the comparisons. However, there are a few comparisons within the 30-39 and 40+ age groups for which the estimated HR for the MI approach is closer to 1 relative to the estimated HR for the SRI approaches. These differences are detected when evaluating the association between hormones and suicide attempts for those without social affirmation and social affirmation for those with hormones but not surgery in the 30-39 year age group and when evaluating the social affirmation for those with hormone or surgery in those 40 years or older. The results are also similar for the MI on the original scale and the log scale. Within each imputation method, we generate 10 sample data sets. Similar point estimates and Table 1. Unadjusted hazard ratios of gender affirmation in binary gender identity group compared among imputation methods. The HRs (IRRs for initial analysis) that are bolded represent those significantly different from one (p<0.05). Reference group = no gender affirmation of any type. Social (No surgery/hormones) = participants had social affirmation but did not have surgery or hormones. Social (Surgery or hormones) = participants had social affirmation and either or both of surgery and hormones.

Imputation Approaches
Initial Analysis slightly narrower 95% confidence intervals were found when we generated 100 data sets (S1 Table).

Discussion
A challenge faced in this analysis of gender affirmation and suicide risk among USTS transgender participants was data missing by design. When evaluating repeated time to event outcomes, it is ideal to use a recurrent survival time model. We were limited in this study due to the fact that the age of each suicide attempt was not collected. To address this challenge, we used three approaches to impute the age at suicide attempt and then fit a recurrent survival analysis with a time-varying covariate for gender affirmation. Similar estimates and conclusions are found for the three imputation methods for the majority of the comparisons (87.5%, 21/24). We observed greater similarity between the two Table 2. Adjusted hazard ratios of gender affirmation in binary gender identity group compared among imputation methods. Model adjusts for age at transgender awareness, race/ethnicity, assigned sex at birth, and education. The HRs (IRRs for initial analysis) that are bolded represent those significantly different from one (p<0.05). Reference group = no gender affirmation of any type. Social (No surgery/hormones) = participants had social affirmation but did not have surgery or hormones. Social (Surgery or hormones) = participants had social affirmation and either or both of surgery and hormones.
Conclusions are consistent with previous survival analysis and analysis of attempt [12], showing negative associations between surgery and suicide attempts in participants who also had social affirmation in a binary gender identity population. But there are some differences between them. For example, in the analysis of attempt rate using zero-inflated negative binomial model, the rate ratios of having surgery versus not in participants who also had social affirmation ranged from 0.7 to 0.9 in participants aged 18-39 and was not statistically significant in participants aged 40 or above (RR 0.9, 0.7-1.0, p = 0.0712). In the recurrent survival analysis using the imputed data, however, the magnitude of association is larger, hazard ratios ranging from 0.4 to 0.6, and statistically significant across all age groups. The cause of the difference is likely to be that the older group was affected the most by lack of temporality in the analysis of attempt rate. This is because, in the zero-inflated negative binomial model, attempts that happened before surgery are lumped together with those that happened after surgery while only the latter is attributable to having surgery.
Another difference from the non-recurrent based analysis is that we additionally found an interaction between social affirmation and hormones which could resolve the contradictory findings on hormones in different age groups. In the zero-inflated negative binomial model [11], we found hormones associated with increased risk in 40+ group but decreased risk in 18-24 age group. However, analysis on imputed data revealed that social affirmation moderates this difference. Specifically, we found hormones without social affirmation associated with increased risk but hormone with social affirmation associated with decreased risk in 18-24 age group. These additional statistically significant findings are probably due to increased statistical power after missing data are imputed and included in the analysis.
We are aware that there are other imputation methods to consider, such as maximum likelihood, and multiple imputation using Cox models. Nevertheless, these alternative methods face some difficulties, including the necessity of specialized software and complex structure of the final analysis model. For instance, Paul Allison has listed several advantages of the multiple likelihood imputation method over MI, including more efficiency, fewer decisions involved, and lack of conflict between the imputation and the analysis model [14]. However, we faced some difficulty implementing the ML method with our data because the PHREG procedure in SAS, which is the procedure we used for all the other analyses, does not impute for missing data. And to our knowledge, there is no ready-to-use software that can impute failure time for recurrent event survival analysis. It is nevertheless an approach worthy of exploration in the future.
We also considered more complex models, such as mixture cure model allowing us to model the probability that participants are cured and will not have any attempts after gender affirmation. The Cox PH model we used assumes that all the recurrent events are identical. This is a commonly used and straightforward model but there are other modeling options that may more accurately reflect the nature of suicide attempts. For example, the accelerated failure time model is a parametric model that assumes the covariates increase or decrease survival time, instead of the hazard [25]. The mixture cure model jointly models a response variable for multiple groups with differential probabilities of susceptibility. This will allow us to assume that participants with certain characteristics or having had certain gender affirmation will be cured and never have any suicide attempts. Our previous use of zero-inflated negative binomial model on count data adopted a similar idea. With the failure time data imputed, all these other modeling options can be explored in the future.
Finally, we recognize that external evidence could inform assumptions for imputation, but such evidence is very limited in this less researched population. Suicide attempt rates by age are available for the general U.S. population [26] but not for transgender population specifically [27]. Yet, suicide patterns may differ greatly between the two groups. It has been reported that the transgender population has a much higher suicide attempt rate than the general population [2]. It is highly likely that the distribution of suicide attempt also differs in these two populations. Therefore, the resultant limitation is that we cannot directly use data in general population to correct for the uncertainty in the suicide probabilities we used for imputation. The USTS survey is the first to provide data on age-specific suicide attempt rates in the transgender population. As more data are gathered, they should be utilized in imputation.
Findings have implications for future collection of suicide attempt data in the transgender population. Optimal suicide attempt measures need to both minimize potential recall bias and capture information of sufficient complexity and completeness to enable the implementation of ideal analytic methods. Future research may consider alternative designs to the cross-sectional survey to potentially reduce the bias. Participant-report measures of suicide attempt history are sensitive questions; thus, care must be taken to balance participant burden and potential upset. Research to validate and compare measures that ask first and last attempt to those which enumerate all attempts in this population will be important for future surveys.

Conclusion
In conclusion, we evaluated three imputation approaches in addressing missing failure times in recurrent-event type of data. In our USTS data of suicide attempts and gender affirmation, parameter estimates are insensitive to imputation approaches. Findings from imputation methods are generally consistent with those from single-event survival analysis and zeroinflated negative binomial model. However, using the imputation approaches, we were able to identify additional associations that were not detected in the zero-inflated negative binomial model and were consistent with clinical expectations. Specifically, we found a statistically significant association between surgery and suicide attempts in 40+ group and an interaction between hormones and social affirmation in 18-24 group, while the zero-inflated negative binomial model failed to do so. These additional findings are likely due to increased statistical power and re-established temporality in the imputed data. Using the imputation approaches, we were able to confirm the findings and identify additional associations that were not detected in the initial analysis using the zero-inflated negative binomial model, highlighting the potential value of using imputation methods in the setting of missing failure times. Findings suggest the need to consider additional analytical approaches in settings with high levels of data missing by design and to conduct further measure validation for future surveys of suicide risk in the transgender population.
Supporting information S1 Fig. Histogram of imputed age by approach (100 imputation datasets). (DOCX) S1 Table. Unadjusted hazard ratios of gender affirmation in binary gender identity group compared among imputation methods with 100 imputation datasets. Reference group = no gender affirmation of any type. Social (No surgery/hormones) = participants had social affirmation but did not have surgery or hormones. Social (Surgery or hormones) = participants had social affirmation and either or both of surgery and hormones. (DOCX) S2 Table. Adjusted hazard ratios of gender affirmation in binary gender identity group compared among imputation methods with 100 imputation datasets. Reference group = no gender affirmation of any type. Social (No surgery/hormones) = participants had social affirmation but did not have surgery or hormones. Social (Surgery or hormones) = participants had social affirmation and either or both of surgery and hormones.